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VP Abstract. This paper is devoted to the validation of the theory of neurogeom- 

f"^ etry of vision, due to Jean Petitot. We focus on the theoretical and numerical 

^^ aspects of integration of a certain hypoelliptic diffusion operator that appears 

in the theory. We provide, at the end, a complete numerical algorithm, fully 

parallellizable. 
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1. Introduction 



^ In his beautiful book [30], Jean Petitot proposes a sub-Riemannian theory for 

the neurogeometry of (human) vision; see also [TU]. For mathematical aspects we 
refer also to our previous paper [B]. The main idea goes back to the 1981 Nobel 
prize of Hiibel an Wiesel, who showed that in the human visual cortex VI, there 
are groups of Neurons that are sensitive to position, and groups of Neurons that 
are sensitive to direction, with connections between them that are activated by the 
image. 

So that it is presumable that VI in fact lifts the images f{x, y) (that are functions 
of two position variables x,y in the plane M^ of the image) to functions over the 
projective tangent bundle PTM.'^. This bundle PTM^ has a contact structure which 
is invariant under the action of the group SE(2) of motions of the plane, that will 
be described precisely just below, and that is denoted by ▲ in vector distribution 
form. 
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2 HYPOELLIPTIC DIFFUSION 

Then, Jean Petitot was naturally led to consider sub-Ricniannian structures over 
PTR- also invariant under the action of SE{2), i.e., motion-invariant Riemannian 
metrics over ▲. 

Since we will be obliged at some point to go to certain stochastic considerations, 
it is convenient to consider the following sub-Riemannian structure in a "control 
system form" 

/ i \ / cos(0) \ / \ 

(1.1) y ]=ul sm{9) ]+v[ = uF + uG, 

where {x,y,9) are the coordinates on PTM.'^; u,v are the controls; F,G are vector 
fields given by the formulae 

(1,2) f.c„w|: + ™ot|, g^A, 



The distribution ▲ is generated by the vector fields F,G defined in (1.2 1. Since 
dim A — 2 and the Hormander conditiorjj holds true at any point, ▲ defines a 
sub-Riemanninan structure on the bundle PTM.'^. The Riemannian metric (•, •) is 
specified via the "controls" u,v, by the fact that F,G form an orthonormal frame 
field. 

The length of an admissible (i.e., everywhere tangent to A) smooth curve F, 



r(t) = {x{t), y{t),e{t)), is therefore L{T) ^ /p y ±2 ^ y2 ^ 02 ^^ ^he sub-Rieman- 
ninan distance between two points of PTMr is the minimum length of admissible 
curves connecting these points. 

In practice, it appears important to introduce a weight parameter a > and the 



corresponding metric -^^(r) = L \J x^ + 2/^ + ck^^ dt. The meaning of the parame- 



ter a will be discussed later (section 2.3.1 1 



Remark 1. 1. The formula (1.1 1 is not global over PTR : this sub-Riemannian 
structure is not trivializable. However, ij we lift it to the group SE(2), which is a 
double covering of PTM.^ , it becomes trivializable and the equation (1.1) becomes 
global. 

2. It was noticed by A. Agrachev that this is the only sub-Riemannian structure 
over PTIR2 -which is invariant under the action of SE{2). 

3. From the theoretical point of view, the weight parameter a is irrelevant: for 
any a > there exists a homothety of the {x,y) -plane that maps geodesies of the 
metric with the weight parameter a to those of the metric with a = I. 

For the same reasons as in Riemannian geometry, minimizing the length L{T) 
is the same as minimizing the energy E{r) — /p(i'^ + y'^ + &^) dt = /p u'^{t)dt -\- 
Jp v'^{t) dt. From the vision point of view, the term /p u^(i) dt is the energy neces- 
sary to activate neurons that are sensitive to position in the plane, while Jp v'^it) dt 
is the energy necessary to activate neurons that are sensitive to direction. 

The idea is that, when reconstructing a corrupted image, the visual cortex 
VI will minimize the energy necessary to activate the neurons that are not activated, 
in the corrupted part of the image. For this, consider an individual interrupted 
level curve f{x,y) = const. It is reasonable to complete it by a geodesic of the 

^The Hormander condition reads: span{F, G, [F, G]} = TqPTR'^ for any q S PTR^. 
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sub-Riemannian structure, that minimizes energy of neurons to be activated. In 
this direction, behind an important theoretical contribution, successful practical 
results were obtained [311 I32j. However, for a corrupted piece of a real image, it is 
not clear how to put in correspondence the non-corrupted parts of the same level 
set. 

For that reason, to reconstruct a corrupted image (the VI inpainting problem), 
it is natural to proceed as follows: in system ( |1.1[ ), excite all possible admissible 
paths in a stochastic way. One gets a stochastic differential equation 

(1.3) 

where ut,Vt are two independent Wiener processes. 
Consider the associated diffusion process: 

dtp 1 2 X-.2 f /r.\ ^ • //i\ ^ \ ^^ 




(1.4) ^ = -A^, A^F^ + G^=^cos(.)-+sin(.)-^^^ , ^^, . 

The operator A is hypoelliptic (satisfies Hormander condition). By the Feynman- 



Kac formula, integrating Equation (1.4) with the corrupted image as the initial 



condition, one expects to reconstruct the most probable missing level 
curves (among admissible). 

To summarize, the analysis of Jean Petitot of the process of reconstruction by 
VI of corrupted images is the following (we refer to [5] for details) : 

• The plane image /(x, y) is lifted to a certain "function" f{x, y, 9) on the 
bundle PTM^. 

• The diffusion process ( |1.4[ ) with the initial condition '/'L q= / is integrated 
on the interval i e [0, T] for some T > 0. 

• The resulting function frp = V'L_7i on the bundle PTM^ is projected down 
to a function fT{x,y), which represents the reconstructed image. 

The lifting procedure should be as follows: the image is assumed to be a smootqj 
function f{x,y). Then, it can be naturally lifted to a surface S in PTM."^ by lifting 
its level curves. At a point {x,y,9) we would hke to set f{x,y,9) = f{x,y) if 
{x,y,9) e S and f{x,y,9) — elsewhere. But this would be nonsense since S has 
zero measure in PTM.'^. Hence in the Petitot model f(x,y) is lifted to a distribution 
f{x,y,9), supported in S, and weighted by f{x,y). 



It turns out that although looking rather simple, Equation (1.4) is not that 
easy to integrate. In particular, the multiscale sub-Riemannian effects are hidden 
inside. We refer to our previous paper T for the computation of the associated 
heat kernel, but the direct use of the kernel appears quickly not very tractable. 
Moreover, the numerical integration starts to be a quite large problem, due to the 
number of points/angles in a reasonable image. To get an idea of the size of the 
full space-discretization of the problem, see Remark |4] below. 

Hence a significant part of the paper is devoted to the numerical integration of 
this equation. 



In fact, it is known by biologists that a certain smoothing procedure is aheady applied at the 
level of the retina [T2ll25ll28| . 
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It turns out that we were led to rather abstract considerations, to get at the end 
a quite simple algorithm, massively parallelizable. This aspect of parallelizability 
seems to fit with the structure of the visual cortex VI, which is suspected, due to 
this structure, to make a lot of parallel computations. 

The paper is organized as follows: in the next section [21 we discuss some prop- 
erties of the group SE{2) and its discrete avatars, together with a discrete version 



of the hypoelliptic diffusion (1.4) 



Section [3] presents our main ideas and shows a first algorithm, which is not 
suitable in our case, mostly for two reasons: 

• The "distributional" nature of lifted images (initial conditions). 

• Compared to our final algorithm, it does not present the capability to treat 
in a different way the corrupted and non-corrupted parts of the image (see 
Section [5]). 

However, this approach could certainly be interesting for more smooth initial con- 
ditions. 

Section [4] presents our final algorithm which, due to the abstract structure of 
the groups under consideration, is massively parallel: It reduces the problem to a 
number of problems of integration of linear differential equations in low dimension, 
completely decoupled. 

Section [5] presents a heuristic improvement of the algorithm in the case where 
we can distinguish between the corrupted and non-corrupted parts of the image. 

Finally, a very short appendix is devoted to basic theoretical concepts and to 
certain technical computations. 

We do not pretend that our results are better than the inpainting algorithms 
existing to day. We just claim that they really seem to validate the theory of Jean 
Petitot and we emphasize the global character of the algorithm. 

We do not discuss the final step of the algorithm (projection) here. There are 
at least two obvious possibilities: projecting by taking either the maximum, or the 
average, over angles. Numerous experiments show that the projection made by the 
maximum provides better results. 

2. The Operators and the groups under consideration 

2.1. Groups. We advise the uninitiated reader to start with our paper [21 and to 
have a look to the appendix at the end of this paper. 

2.1.1. The group of Motions SE(2). The group law over the Lie group SE{2) has 
multiphcation law (^2,62) • (-'^1, ^1) — {^2 + Re2^iiOi + ^2), where 

cos(6') sin(6l) 



^ ~\y )' ^^^ \-sm{e) cos{9) 

Re is the rotation of angle 9 in the X-plane. 

The (strongly continuous) unitary irreducible representations of SE(2) are well 
known. However for analogy with the next section, we need to recall basic facts. 
For a survey, see [Ml [33] . 

As representations of a semi-direct product, they can be obtained by using 
Mackey's imprimitivity theorem, and therefore, they are parametrized by the orbits 
of the (contragredient) action of rotations on the X-plane, i.e., they are parametrized 
by the half lines passing through the origin. Additionally, corresponding to the ori- 
gin, there are the characters of the rotation group S^ that do not count in the 
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support of the Planchcrel's measure. Finally, the representation x\ corresponding 
to the circle of radius A acts over L'^{S^), and is given by 

where v\ is a vector of length A in M'^ and (•, •) denotes the standard Euclidean 
scalar product over M? . 

Note that SE{2) is far from being maximally almost periodic, since all its finite 
dimensional unitary irreducible representations are given by the characters of S^ 
only. 

Let us also recall [21 [6l [131 [14] that our sub-Riemannian heat kernel corresponding 



to the hypoelliptic Laplacian ^A defined in (1.4) is given by 



(2.1) P,iX, e)^\ j {^ e^'l^ce,, (^9, ^) , xx{X, 9) ce„ (^0, ^ 

^ e^-'/se Je, M , xa(X, ^) seje, ^"j \ AdA, 

where the functions se„ and ce„ are the 27r-periodic Mathieu sines and cosines, and 
a^j = — ^j^ — a„(^), b^^ — —^ — &„(^), with a„,6„, the characteristic values for 
the Mathieu equation. 



Then the solution of the diffusion equation ([1.4 ) with the initial condition V'L-n^ 



'00 is given by the right-convolution formula: 

(2.2) 

i;{t,X,9)^e'^MX,9)^MX,9)*Pt{X,9)^ J Mg)Pt{g-^ ■ (X,9)) dg, 

SE{2) 

where g G SE{2) and dg is the Haar measure: dg — dxdyd9. Unfortunately, in 



practice formula (2.2 1 is not very tractable, for several reasons (in particular, due 
to the slow convergence and the difficulty to find good computer implementations 
of Mathieu functions). 

2.1.2. The group of discrete motions SE{2, N). In the next section, a discrete-angle 
version of the above diffusion equation will appear naturally. It corresponds to the 
group of motions SE{2) restricted to rotations with angles ^. This group is 
denoted by SE{2, N). It has very special features: it is maximally almost periodic, 
and all its unitary irreducible representations are finite dimensional (it is a Moore 
group, see |5D] for details), although it is not compact. This follows in particular 
from [TBJ 16.5.3, page 304]: it is the semi-direct product of a compact subgroup K — 
Z/iVZ, and a normal subgroup V isomorphic to M^, each element of V commuting 
with the connected component of the identity (which in this case is V itself). To 
simplify, we set Rk = R2ki, . 

N 

Besides the characters of K^ the representations coming in the support of the 
Plancherel measure are parametrized by (A, v) G IR+ x [0, ^[~SE{2, N), the duaP 
of SE{2, N). They act on C^ and are given by 

(2.3) XxAX,r) = diagfc(e^<^^-«^^>)^^ 



The natural "dual topology" of SE{2, N) is that of a cone, that consists of considering 
' X [0, ^] and identifying (A,0) with (A, ^). 
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where diag^, (^gnVA.,^.fifcX> j jg ^^j^^ diagonal matrix with diagonal elements e'^^-^'^-^'''^^ , 
V\,i, = (Acos(i/), Asin(i/)), k — 1,. . . ,N, and S is the shift matrix over C^ (i.e., 
Sck = Ek+i ioT k — 1, . . . , N — 1, and Scn = ei). 

2.2. The semi-discrete diffusion operator. 

2.2.1. Semi-discrete versus continuous. Firstly, we show that a certain semi-discrete 
(discretization with respect to the angle) model of the diffusion is compatible with 
the limit continuous model. For the considerations in this section, one may refer to 
the paper dj. 



The diffusion equation ( 1.4 ) comes from the stochastic differential equation ( 1.3 1, 



where ut, Vt are two independent standard Wiener processes. The diffusion equation 



(1.4) is the associated Fokker-Planck (or Kolmogorov forward) equation for (1.3). 

From the image processing point of view, integrating the diffusion is equivalent 
to excite all possible admissible paths, in a stochastic way. In fact, in the real 
structure of the VI cortex, not all angles are represented, but a small 
number only. This number is denoted by N. 

Therefore, it is natural to consider the following stochastic process with jumps 



evidently connected with the stochastic equation (1.3): 



in which 6t is a jump process. Set A^v — {K.j), h j — 0, . . . , N — 1, where 



P [9t = ej\9o ^ e,] . . , . , y^ 



A,; J = lim for i ^ j, X^j = - 2^ A^, 

The matrix A^r is the infinitesimal generator of the jump process. 

We assume Markov processes, where the law of the first jump time is exponential, 
with parameter /3 > (that will be specified later on). The jump has probability 
i on both sides. Then we get a Poisson process, and the probability of k jumps 
between and t is 

P [/fc jumps] = ^e-'5*. 

So that: 

P [0t = e.±i|0o = e,] = l{Pt + kt" + o{t'))e'^\ 

P [dt = e,±r.% = e.] = 0(t")e-'^*, n = 2, 3, . . . , 

with the convention that e^ is modulo N. 

Hence Ai^i±i = ^/3 and Xi^i = — /3. All other elements of the matrix Ajv are equal 
to zero. Then, the infinitesimal generator of the semi-group associated with the 
stochastic process {zt,Ot) is of the form: 

(2.5) Cn^{z, a) = {A^),{z) + {Kn^{z, e,)),, 

where ^j(z) = ^(z, e^), z ~ (x,y), and 

If , . d ... 9^^ 



(2.6) (A*),(z) = A*(z,e,) = -(cos(e,)— -Hsin(e,)— ) *(z,e, 



HYPOELLIPTIC DIFFUSION 



(2.7) 



(AA,*(z,e,)), 



n-l 



P, 



^ A.,,M/,(z) = |(*,_i(z) - 2*,(z) + vI/,+i(z)). 



JV\2 



Therefore, if we set /3 = (5^)", we get 
l*,_i(z)-2*,(^) 



(Ajv*(z,e,)), 



^,+i(z) 1 5 



'27r\2 



2^,,*(^,e,0 + O^^ 



At the hmit N ^ 00, from formulae (2.5) ~ (2.7) we get the second order differential 
operator 

2 02 



C^iz, 



cosW-+sm(^)- 



502 



1 



^{z,e) = -A^{z 



that is, the operator of our diffusion process (1.4). However the exact Fokker-Planck 
equatiory with number of angles N ^ 00 contains the parameter /3: 



(2.8) f:(M) 



-f cos(er)^ +sin(er)— ) Vr(i, ^)+ 



/3 



{A-iit, z) - 2Mt, z) + i^r+iit, z)), r = 0, . . . , iV - 1. 



2.2.2. The semi-discrete heat kernel via the GFT. The GFT (generalized non- 
commutative Fourier transform, see [51 and Appendix here) transforms our hypoel- 
liptic diffusion equation into a continuous sum of diffusions with elliptic right-hand 
term, and the summation over SE{2, N) is with respect to the Plancherel measure 
\dXdv. We compute the semi-discrete heat kernel via the GFT, just as in '^, and 
we get a similar but simpler formula than for the "continuous" heat kerne (2.1) in 
the case of the group SE[2): 

Formulae (6.1) and (6.2) for the direct and the inverse GFT show that, if we set 

(2.9) Ia,^ = Aat -diagfc(A2cos2(efc - v)), 

we get the following expression for the "semi-discrete" (or "jump") heat kernel: 

(2.10) Dt{z,er)= I trace('e^^-*-diagJe'<^^.-^^'=^>)s"')AdAdz^. 

5B(27/V) 

Note also that it is a closed formula similar to the one in the Heisenberg case 
(explicit usual functions modulo a Fourier transform) . These are the only cases we 
know for noncompact groups, where the kernel is obtained with such a formula. 

2.2.3. A direct way to compute the semi-discrete heat kernel. We start from our 
semi-discrete heat equation (2.8). The heat kernel is the convolution kernel associ- 
ated with Dt{z, Cr), the fundamental solution of Equation (2.8). 

Applying the ordinary Fourier Transform with respect to the space variable z 
to ( |2.8[ ), we get the ordinary linear differential equation 

(2.11) ^5t{z) ^ Ax,Mz), 

dt 



Note that here, due to self adjointness, the Kolmogorov-backward equation is the same as the 
Kolmogorov- forward one (Fokker-Planck). 
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where Dt{z) ~ {Dt{z, ei), . . . , Dt{z, cn)), with the initial condition obtained from 
the Dirac delta function at the identity by the ordinary Fourier Transform: 

5o(z)=<5jv = (0,...,0,l). 

Then, taking the inverse ordinary Fourier Transform with respect to the space 
variable, we get a second expression for the heat kernel: 

(2.12) Dt{z,er)^ l(e^^-"'6N\ e^<^^-">AdAd:/, 

where the subscript r means the r-th coordinate of a vector. 



Remark 2. Looking at the formulae (2.10) and (2.12), it is not clear that these ex- 



pressions are identical. The proof of this fact is given in the appendix (section 6.3) 



Although less direct, we prefer formula (2.101, that reflects more the structure of 

SE{2,N). 

2.3. Two important points. 

2.3.1. The weighting of the metric. Consider the diffusion process 

(2.13) f = ^A.^, A„^F^+oG^^(cos(.)A+3in(.)|)%«^ 



with the weighting coefhcient a > 0. The hypoelliptic operator A defined in ( |1.4[ ) is 
the case a = 1. Although the coefficient a is theoretically irrelevant (see Remark [l] 
item 3), it plays an essential role in practice. 

On the same way, for the semi-discrete operator, we have 

AWV^(i,^) = 

cos(ei)^ +sin(ei)— j ipi{t, z) + f3{;tpi^i{t,z) - 2ijji{t, z) + i/^i+iit, z)) = 

cos(e.)^+sm(eO-j + p[-) ^j^.(M) + o(-) 



as A'^ — >■ 00. Comparing the above formula with (2.8 1, we have the relation: 

(2.14) "-Klv)- 

It appeared clearly in all our experiments that TV = 30 is always enough (it brings 

nothing visible in the experiments to take N > 30). 

Both parameters N and /3 have a physiological meaning. Therefore, 
the coefficient a of the limit behavior [N -^ od) can certainly be obtained from 
physiological considerations. 

2.3.2. The effect of the projectivisation on the kernel. This is related to the known 
phenomenon called "chirality of the pinwheels" : the pinwheels are cells in some 
pictures that are obtained by observing cuts of the VI cortex. They put in evidence 
the sensitivity of certain neurons to directions, and they show clearly that certain 



groups of neurons have spatial chirality. This could be reflected in formula (2.15) 
below. 

As we said, the neurons are not sensitive to the angles themselves, but to di- 
rections only (i.e., the angles are modulo tt and not 27r). It is the reason why we 
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work in the projectivisation PTR^ of the tangent bundle of M^, in place of SE{2) 
itself. From the discrete point of view, it means that if N is the number of values 
of directions (not angles) , the angle-step is in fact -^ . 

Also, as it was explained in [6^, the sub-Riemannian structure over PTM.'^ itself 
is not trivializable. This is not a problem for the operators A and A'^', since the 
functions cos(6')^, am{9)^ and sm{29) are vr-periodic. Hence, details relative to 
this projectivisation are omitted in the following sections. 

From the point of view of the heat kernels, in fact, if pt denotes the heat kernel 
over PTM.'^ (which is not a group convolution kernel anymore, since PTM^ is not a 
group), we have the formula 

(2.15) pt{{x, y, 9), {x, y, 9)) = Pt{{x, y, 9)-' ■ (x, y, 9))+Pti{x, y, 9)-' ■ (x, y, + 7r)), 

where the inverses and products are intended in the group SE{2). 



A similar formula holds for the semi-discrete kernels dt and Dt (from (2.10l) for 
even N . 

3. A PRE-ALGORITHM 



As we said in section 2.1.2 the group SE(2,N) is maximally almost periodic. 
It follows from the expression ( |2.3| ) of the unitary irreducible representations that 
the Bohr-almost periodic functions f{x, y, r) are just those such that the functions 
fr{x,y) = f(x,y,r), r = 1, ... ,7V, are Bohr-almost periodic over M? in the usual 
sense. We call AP[N) the set of almost periodic functions on SE{2,N), and we 
identify the elements of AP{N) to C^-valued functions whose components are 
almost periodic over M^, i.e., functions that are uniform limits of trigonometric 
polynomials in the two variables {x,y). 

These functions are dense among continuous functions over any compact subset 
of SE{2, N), and AP{N) is a good candidate for the space of solutions of our heat 
equation: exactly as for the usual heat equation (see [TTl pp. 144-146] for instance), 
the uniformly bounded solutions of our heat equation with initial conditions in 
AP(N) remain almost periodic over SE{2, N) uniformly in time. 

Functions coming from practical images have support in a bounded subset of 
SE{2,N). They can be (after smoothing and lifting) approximated uniformly over 
this bounded subset by trigonometric C^-valued polynomials Q{x,y) with compo- 
nents Qr{x,y): 

k 
(3.1) Qr{x,y) = ^a,,A.,^,e^(^'^+^'^^ 

1=0 

The vector space of such C^-valued polynomials, for a fixed finite number of distinct 
values of Xi,fJ.i, i.e., u — {Xi,fii) £ K, a, fixed finite subset of M^, is denoted by 

SE{2,N,K). 



A trivial computation shows that semi-discrete hypoelliptic equation (2.8) re- 
stricts to SE{2,N,K). It becomes 

(3.2) ^^-l (^' cos(0,) + M. sin(0,)) 4 + f (<+i " Saj, + al_,] , 

or, equivalently, 

dA^ 1 2 

(3.3) — = --diag^:{X,cos{9k) + ^i^sm{9k)) A' + AnA\ 
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where A^ is the vector {a\, . . . , a^y) and the matrix Ajv is the infinitesimal generator 



of the jump process with parameter /?. The system of differential equations (3.2 1 is 



equipped with the initial condition a^(0) — ar.\-^^- with ar,x.^^. from (3.1). 
The following theorem holds. 

Theorem 1. For any almost periodic polynomial initial condition in SE{2, N, K), 
the integration of the diffusion equation reduces to solving a finite set of indepen- 
dent linear ordinary differential equations of dimension N. Any continuous initial 
condition over a compact subset of SE{2, N) can be uniformly approximated by such 
a polynomial. 



Differential equations (3.2 1 over C^ are not hard to tract numerically. They 
have an elliptic right hand term, and the Crank-Nicolson method (discussed in the 
next section) is recommended. 

In fact, SE{2, N, K) is the space of a unitary representation of SE{2, N) which 
is not irreducible but splits into the direct sum of irreducible representations: 



SEi2,N,K)^^SE{2,N,{Lo}). 



ujeK 



This fact suggests that we can use our knowledge of the dual SE{2, N) to reduce the 
computations. Actually, it is easy to check that for lj^uj £ K such that uj = RrUJ, 



if we call A and A the corresponding solutions of Equation (3.3 1 in SE(2, N, {uj}) 
and SE{2,N,{uj}), respectively, we get: 

— = --diagj.(AjCos(6'fc) + AijSin(6'fe)) A + A^A, 

1 Qf A 1 

— ^ = --diagfc(A, cos(6ifc) + fn sin(6ife))^5M + AnS^'A. 

So that it is not necessary to compute all the resolvents relative to each unitary 
irreducible representation. It is enough to do it for those corresponding to cj G 

SE{2rN). 

Theorem 2. In theorem^ if some points of K can be deduced one from the other by 
elementary rotations Rr then it is enough to compute the resolvents corresponding 
to a single among these points. 

This method could be very efficient in a general setting. The last considerations 
show that it is of numerical interest to put in K as much as possible points in the 
same orbits under the elementary rotations. 

In our vision problem, the initial conditions are certainly very far from being 
almost periodic. Hence we preferred a more direct method, which is however closely 
related with this one. This is explained in the next section. 

4. Final algorithm 

In fact, images being given under the guise of a square table of real values (the 
grey levels), we have chosen to deal with periodic images over a basic square, and 
to discretize with respect to the naturally discrete (x, y) variables. We take a mesh 
of Mj. X My points on the (x, 2/)-plane, and the number of angles is N . 
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Remark 3. In practice, in all the results we show, M^ = My = M = 256, and the 
number of angles N w 30. We don't see any significant improvement for 
larger N. 

Therefore, a function ip{x,y,0) over SE{2,N) is approximated by a M^ x N 
table (V'f i), P = 1, ■ ■ • 7 -^7 k,l — 1, . . . , A/. Hence, without loss of generality, we set 
for the discretization steps Ax ~ Ay = ^/M and the mesh points are Xk — ^^, 
yi = -jf7- In this section, due to the periodicity, the upper index takes natural 
values modulo iV and the lower indices take natural values modulo M. 

Remind that the semi-discrete diffusion equation over SE{2,N) is (2.8). For 

^T JL 

dx ' dy 



numerical solution of this equation we replace the differential operators -^, ^^ by 
the discrete operators D^, Dy, which act on C*^ 



/,P A - "^fc+l.' "^fe-l-' - ^^/,P _ .l,P 



^xKVk.l) — ^ _^ ^9 \Vk+l,l 'f^k-l.U' 

Xk+1 Xk—l ^ 

r^ ( IV \ '^k,l+l~'^k.l-l y^ / IP IP \ 

Then we get the diffusion equation with totally discretized (i.e., with discretized 
space and angle variables) operator D, which acts on C^ (g) C*^ ~ (J^wm . 

dTpl, 1 
(4.1) ^=2^«;)' ^here 

Di^k,) = {cos{e,)D, + ,\n{e,)Dy)\rk,) + /^IV-^^' " '^^h + <f)- 



The initial condition for ( |4.1[ ) is the discrete analog of the function f{x,y,9) ob- 
tained after the lift of the original image f{x,y). 



Remark 4. For M = 256, TV = 30 ( |4.1[ ) is a fully coupled linear differential equa- 
tion in M^, K — 1,996,080. Applying an implicit or semi-implicit finite difference 
scheme, one need to solve a system of K^ ~ 3.6 x 10^^ linear algebraic equations. 



As previously, it is natural to apply the Fourier transform over the abelian group 
(Z/AfZ) , which can be computed exactly by the standard FFT (Fast Fourier 
Transform) algorithm. Let us denote by ip the Fourier transform of ip: 

^ 1 *^ / 



'ik-l)ir-l) , il-l)is 



r,s — 1 ^ ^ 

A straightforward computation shows that 

OM.^ - ^^Msin^2n^yjl, i?;(^0 = ^yMsin(^2.^) V^^,, 

and the Fourier transform maps the operator cos{ep)D,j. -\- sm{ep)Dy to the multi- 
plication operator -0^; — > i^/M a^^, li^k i' where 



a^i = cos(ep) sinf 2-7r J + sin(ep) sinf 27r— — : 
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Hence, the diffusion equation (4.f I is mapped to the following completely decoupled 
system of M^ linear differential equations over C^: 

(4.2) ^=A,,,^fc,,, ^,,; = (^i,,,...,V^^,), 

at 

where Ak.i = ^{-^n ~ ^^diagp(a^ ;)^), and A^ is an almost tridiagonal N x N 
matrix: due to periodicity in p, An contains two extra non-zero elements in the 
right-up and left-low corners. 



Therefore the solution of (4.2) is 
(4.3) V^fe,i(i) = e*^'=.'^fc/0), 



where the initial function ^pk,i{0) is known. Finally, the solution of (4.11 is the 
inverse Fourier transform over (Z/AfZ) of tp^ ^(i), obtained using the inverse FFT. 



Theorem 3. The solution of Equation (4.1) can be computed exactly by solving 



in parallel KP linear differential equations in dimension N . 

For instance, for M = 256 and N ~ 30, this is solving 65536 linear differential 
equations in dimension 30. 

Remark 5. 1. The complexity of the FFT used twice in our algorithm (includ- 
ing the inverse transform) is negligible in front of the complexity of the numerical 
integration of the decoupled linear differential equations. 



2. The discretized diffusion (4.1) can be numerically integrated with the semi- 
explicit Crank- Nicols on scheme. The convergence of this scheme for the considered 
class of equations and some estimations are well known; see e.g. (HI chapter 5] . 
Note that the matrices A^^i are tridiagonal plus two terms in the right-up and left- 
low corners (coming from An due to periodicity). An effective algorithm for solving 
such linear system is suggested in 1 . 



3. Formula (4.3) for solutions of the linear differential equations (4.2) has an 
obvious advantage. Once M, N and a (or, equivalently, /3) are fixed, the matrices 
Ak^i are universal, and therefore their exponentials can be computed once for all. 
Moreover, using a time step r and the formula e""^ *=' = (e"^ *='') , it is enough to 
compute only the exponential e'^ ''•' . 

4. In fact, it is not necessary to compute M"^ exponentials of matrices: it is 
enough to compute them at the points k,l that fall in the dual space SE{2,N). 
This is much less. For large N, is number is of order M in place of M^ . Here, the 



structure of SE{2, N) is crucial, and specially formula (2.10) 



5. Preservation of non-corrupted parts 

It has to be noticed that: 

• The treatment of images up to now is essentially global. 

• It is not necessary to know where the image is corrupted. 

Presumably, the visual cortex VI is also able in some cases to detect that the 
image is corrupted at some place, and to take this a-priori knowledge into account. 
We tried to investigate a method that would improve on our algorithm in this 
direction. We suggest (essentially heuristically) the following procedure. 



HYPOELLIPTIC DIFFUSION 



13 





Figure 1. The initial corrupted image (left) and the image re- 
constructed via diffusion (right). 
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Figure 2. The initial corrupted image (left) and the image re- 
constructed via diffusion (right). 



5.1. Motivation. One natural idea would be to iterate the following steps: 

• Lift the plane image to PTM.'^. 

• Integrate the diffusion equation for some fixed small r. 

• Project the solution down to the plane. 

• Restore the non-corrupted part of the initial image. 

In practice, this idea does not work at all. The main reason is that the diffusion 
acts on both the corrupted and the non-corrupted parts, and there is not good 
coincidence at the frontier after the restoration. 

All variations of the previous idea, at the level of the plane image do not provide 
acceptable results. From what we conclude that it is necessary to proceed the 
restoration of the the non-corrupted part at the level of the lifted image, i.e., on 
the bundle PTW^. 
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5.2. Description of the "preservation" procedure. Assume that points (x, y) 
of the image under reconstruction are separated into the set G of "good" (non- 
corrupted) and the set B of "bad" (corrupted) points; /(x, y) = for all [x, y) G B 
and f{x,y) > for all {x,y) e G. We would like to preserve points of the set G 
from the effect of diffusion, while on the set B the diffusion should be proceeded. 

The general idea of the "preservation" procedure is to "mix" the solution ipix, y, 9, t) 
of the diffusion equation with the initial function ip[x,y,9,Q) — f{x,y,6) at each 
point {x, y) G G. The "mixing" is fulfilled many times during the integration of our 
equation. 

Namely, split the segment [0, T] into n small intervals with the mesh points 
ti — IT, T — T/n, i — 0,1, ... ,n, and proceed the integration of our equation on 
each [ti,ti^i\ with the initial condition 

I j ■ip{x,y,0,U), a {x,y) e B, 

'*=*' \^a{x,y,ti)i){x,y,9,ti), a {x,y)&G, 

where the function 4'i^j J/? ^i U) is computed after integration on the previous inter- 
val, and the factor a{x,y,ti) is defined by 

e/i(a;,y,0) + {I - e)h{x,y,ti\ 



a{x,y,U) 



h{x, y, t) — max V'(a;, y, 9, t), 



h{x,y,U) 

where < e < 1. After that we obtain the function ip{x,y,9,ti^i) and repeat the 
procedure on the next interval. 

Note that this procedure essentially depends on two parameters that can be 
chosen experimentally: the natural number n (number of treatments) and the co- 
efficient e, which defines the strength of each treatment. 



ti-' -^i 

IB,.. ^ 
■ ■■9D- 



■■■■Kai 
iHBaaHi 



Ik .::^l 



tW JBUNBaiBB 
■ l JBBBnHBB 

IS.' j»BB'/HBB 
■( iBBBBCaBBB 
IBy.:44BBBBBBBB 
IK:"::S»:^BBBBBBBB 
■L^aeaBBBBBBB 




Figure 3. The initial corrupted image (left, 32% damaged) and 
the image reconstructed via diffusion with the "preservation" pro- 
cedure (right). 
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Figure 4. The initial corrupted image (left, 37% damaged) and 
the image reconstructed via diffusion with the "preservation" pro- 
cedure (right). 
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Figure 5. The initial corrupted image (left, 37% damaged) and 
the image reconstructed via diffusion with the "preservation" pro- 
cedure (right). 

6. Appendix 

In this appendix, we collect shortly a few results needed for the understanding 
of the paper. 

6.1. The Generalized Fourier Transform (GFT). Given a locally compact 
unimodular topological group G of Type FJ the dual G is the set of (strongly) 
continuous complex unitary irreducible representations (xg, Hg) of G. For g ^ G, 
Xg- G ^ U{Hg), the unitary group of the complex Hilbert space Hg. The space 
of complex L^ functions over G with respect to the Haar measure is denoted by 



We do not say what type I means. We just need here to know that both SE{2) and SE(2, N) 
are type I. For instance, any connected semi-simple or nilpotent group is type I. 
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L^{g, dg). The generalized Fourier transform (GFT) of / £ L'^{g, dg) is defined bjj^ 
(6.1) m^ f f{g)Xsi9-')dg. 



The operator f(g) is Hilbert-Schmidt over Hg and there is a measure dg over 
G, called the Plancherel measure, such that the GFT is an isometry to L^{g,dg), 
where the space L^ (g, dg) is the continuous Hilbcrt sum of the spaces of Hilbert- 
Schmidt operators over the spaces Hg, with respect to Plancherel's measure. As a 
consequence, we have the inversion formula: 

(6.2) f{g)= Imx-Mdg- 



The GFT is a natural extension of the ordinary Fourier transform over abelian 
groups and it has all the corresponding properties, such as: mapping convolution 
to product, etc. (see [5]). 

As in the case of the usual heat equation over M", we use it to solve our heat 
equations over the groups SE{2) and SE{2,N). 

6.2. Bohr Compactification and Almost Periodic Functions. The Bohr com- 
pactification of a topological group G is the universal object {G^ ,a) in the category 
of diagrams: 

a-.G-^H, 

where ct is a continuous homomorphism from G to a compact group H. If the 
mapping ct: G — >■ G'' is injective, G is called maximally almost periodic (MAP). 

The set AP{G) of almost periodic functions over G is the pull back by a of the 
set of continuous functions over G*" . 

The group G is MAP iff the continuous unitary finite dimensional representations 
of G separate the points. A connected locally compact group is MAP iff it is the 
direct product of a compact group by M". For instance, SE{2,N) is MAP, while 
SE{2) is not. 

A continuous function / e AP{G) iff its right (or left) translated form a relatively 
compact subset of E{G), the set of bounded continuous functions over G, iff it is 
a uniform limit of coefficients of unitary irreducible representations of G. If G is 
MAP, AP{G) is dense in the space of continuous functions over G, in the topology 
of uniform convergence over compact sets. 

For duality over MAP groups, see [H] and the book [50]. For introduction to 
almost periodic functions, see the original paper [35] and the nice exposition [16; . 



6.3. Proof that expressions (2.10) and (2.121 are identical. In this section. 



all integer indices take values between 1 and N, and addition is always modulo N. 



As usual, the integral below is well defined for / G L^{G,dg) only, but is extended by 
continuity. 
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Set Mx,j{t) — e^ >-'"*■ , where the matrix Ax,j is defined in (2.9). Then we need 
to estabhsh the identity of the two expressions: 

Dl{z,er)= I trace(AfA,^(t)-diags,('e'<^^.--"*^>)s"')AdAdi/, 

S£(2JV) 

A'(2, e.) =y"(MA,,(t)<5jv), e'<''^-""> AdAdi.. 

R2 

Set 

Ml^uit) = e^^'-*, Al^u - Aat - diagfe(A2cos2(efc+, - v)). 
The foUowing fact is cruciah 

(6.3) S-'^MxAt)S^' = Ml,{t), 

where S is the shift matrix over C^ (i-e., Sck — Ck+i for fc = 1, . . . , A'^ — 1, and 



Scn — ei). The equahty (6.3) follows from S ^A\ ^5*" — Ax^v, a consequence of 



5 ^'Aj^S^ — An and of the general relation 

(6.4) [b DO )n,m ^ t>n—r,7n—ri 

which holds true for any N x N matrix B. 
An immediate computation shows that 



.5) D\{z,er)^ I ^(AfA,.(t)-diagfe(e'<^--^'=^>))^^^^AdAdz/- 

.— - — -. n ' 

SE{2,N) 



SE{2,N) 



On the other hand, using the relations (6.3) and (6.4), we have: 



SE{2,N) 



DUz,er) = l{MxAt)SN)^e'<''>-"''^XdXdjy = /(MA,.(i))^_^ e^<^^"-^>AdAdi. = 
/ E(A/A%(i)).,^.e*<^^---^^>AdAd^ = 

J-^ n 

;(2,7V) 

/ ^(5-"AfA,.(t)5")^_^e»<^^-«-"^>AdAdi. = 

SE(2,N) 

I 5:(A/,,.(t)),,_„^^_„e^<^--^--)AdAd.. 

SE(2,N) 

Finally, changing n for —n and observing that N -\- n = n, we get the following 
expression: 

Dl{z,er)= f ^(Af,,.(t)),^^^^„e^<^^.-«"^)AdAd^, 

SE{2,N) 



which is exactly (6.5) 
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